Getting started

The packages used in this tutorial can be installed with the following commands:

install.packages("tidyverse")
install.packages("sf")
install.packages("raster")
install.packages("spData")
install.packages("devtools")
devtools::install_github("Nowosad/spDataLarge")

Raster-vector interactions

This section focuses on interactions between raster and vector geographic data models. It includes four main techniques:

  • raster cropping and masking using vector objects;
  • extracting raster values using different types of vector data;
  • raster-vector conversion.

The above concepts are demonstrated using data used in previous chapters to understand their potential real-world applications.

Let’s start by loading the packages we will require for this section of the workshop:

library("tidyverse")
library("sf")
library("raster")
library("spData")
library("spDataLarge")

Raster cropping

Many geographic data projects involve integrating data from many different sources, such as remote sensing images (rasters) and administrative boundaries (vectors). Often the extent of input raster datasets is larger than the area of interest. In this case raster cropping and masking are useful for unifying the spatial extent of input data. Both operations reduce object memory use and associated computational resources for subsequent analysis steps, and may be a necessary preprocessing step before creating attractive maps involving raster data.

We will use two objects to illustrate raster cropping:

  • A raster object srtm representing elevation (meters above sea level) in south-western Utah.
  • A vector (sf) object zion representing Zion National Park.

Both target and cropping objects must have the same projection. The following code chunk therefore not only loads the datasets, from the spDataLarge package, it also reprojects zion:

srtm = raster(system.file("raster/srtm.tif", package = "spDataLarge"))
zion = st_read(system.file("vector/zion.gpkg", package = "spDataLarge"))
zion = st_transform(zion, projection(srtm))

We will use crop() from the raster package to crop the srtm raster. crop() reduces the rectangular extent of the object passed to its first argument based on the extent of the object passed to its second argument, as demonstrated in the command below (note the smaller extent of the raster background in panel (B) below):

srtm_cropped = crop(srtm, zion)

Related to crop() is the raster function mask(), which sets values outside of the bounds of the object passed to its second argument to NA. The following command therefore masks every cell outside of the Zion National Park boundaries:

srtm_masked = mask(srtm, zion)

Changing the settings of mask() yields different results. Setting updatevalue = 0, for example, will set all pixels outside the national park to 0. Setting inverse = TRUE will mask everything inside the bounds of the park (see ?mask for details) (panel (D)).

srtm_inv_masked = mask(srtm, zion, inverse = TRUE)
Illustration of raster cropping and raster masking.

Illustration of raster cropping and raster masking.

Raster extraction

Raster extraction is the process of identifying and returning the values associated with a ‘target’ raster at specific locations, based on a (typically vector) geographic ‘selector’ object. The results depend on the type of selector used (points, lines or polygons) and arguments passed to the raster::extract() function, which we use to demonstrate raster extraction. The reverse of raster extraction — assigning raster cell values based on vector objects — is rasterization, described in Section ‘Rasterization’.

The simplest example is extracting the value of a raster cell at specific points. For this purpose, we will use zion_points, which contain a sample of 30 locations within the Zion National Park. The following command extracts elevation values from srtm and assigns the resulting vector to a new column (elevation) in the zion_points dataset:

data("zion_points", package = "spDataLarge")
zion_points$elevation = raster::extract(srtm, zion_points)
Locations of points used for raster extraction.

Locations of points used for raster extraction.

Raster extraction also works with polygons. Polygons tend to return many raster values per polygon. This is demonstrated in the command below, which results in a data frame with column names ID (the row number of the polygon) and srtm (associated elevation values):

zion_srtm_values = raster::extract(x = srtm, y = zion, df = TRUE)

Such results can be used to generate summary statistics for raster values per polygon, for example to characterize a single region or to compare many regions. The generation of summary statistics is demonstrated in the code below, which creates the object zion_srtm_df containing summary statistics for elevation values in Zion National Park:

group_by(zion_srtm_values, ID) %>% 
  summarize_at(vars(srtm), list(~min, ~mean, ~max))
## # A tibble: 1 x 4
##      ID   min  mean   max
##   <dbl> <dbl> <dbl> <dbl>
## 1     1  1122 1818.  2661

The preceding code chunk used the tidyverse to provide summary statistics for cell values per polygon ID. The results provide useful summaries, for example that the maximum height in the park is around 2,661 meters (other summary statistics, such as standard deviation, can also be calculated in this way). Because there is only one polygon in the example a data frame with a single row is returned; however, the method works when multiple selector polygons are used.

The same approach works for counting occurrences of categorical raster values within polygons. This is illustrated with a land cover dataset (nlcd) from the spDataLarge package, and demonstrated in the code below:

zion_nlcd = raster::extract(nlcd, zion, df = TRUE, factors = TRUE) 
dplyr::select(zion_nlcd, ID, levels) %>% 
  tidyr::gather(key, value, -ID) %>%
  group_by(ID, key, value) %>%
  tally() %>% 
  tidyr::spread(value, n, fill = 0)
## # A tibble: 1 x 9
## # Groups:   ID, key [1]
##      ID key   Barren Cultivated Developed Forest Herbaceous Shrubland
##   <dbl> <chr>  <dbl>      <dbl>     <dbl>  <dbl>      <dbl>     <dbl>
## 1     1 leve~  98285         62      4205 298299        235    203701
## # ... with 1 more variable: Wetlands <dbl>
Area used for continuous (left) and categorical (right) raster extraction.

Area used for continuous (left) and categorical (right) raster extraction.

Rasterization

Rasterization is the conversion of vector objects into their representation in raster objects. Usually, the output raster is used for quantitative analysis (e.g., analysis of terrain) or modeling. As we saw earlier, the raster data model has some characteristics that make it conducive to certain methods. Furthermore, the process of rasterization can help simplify datasets because the resulting values all have the same spatial resolution: rasterization can be seen as a special type of geographic data aggregation.

The raster package contains the function rasterize() for doing this work. Its first two arguments are, x, vector object to be rasterized and, y, a ‘template raster’ object defining the extent, resolution and CRS of the output. The geographic resolution of the input raster has a major impact on the results: if it is too low (cell size is too large), the result may miss the full geographic variability of the vector data; if it is too high, computational times may be excessive. There are no simple rules to follow when deciding an appropriate geographic resolution, which is heavily dependent on the intended use of the results. Often the target resolution is imposed on the user, for example when the output of rasterization needs to be aligned to the existing raster.

To demonstrate rasterization in action, we will use a template raster that has the same extent and CRS as the input vector data cycle_hire_osm_projected and spatial resolution of 1000 meters:

cycle_hire_osm_projected = st_transform(cycle_hire_osm, 27700)
raster_template = raster(extent(cycle_hire_osm_projected), resolution = 1000,
                         crs = st_crs(cycle_hire_osm_projected)$proj4string)

Rasterization is a very flexible operation: the results depend not only on the nature of the template raster, but also on the type of input vector (e.g., points, polygons) and a variety of arguments taken by the rasterize() function.

To illustrate this flexibility we will try three different approaches to rasterization. First, we create a raster representing the presence or absence of cycle hire points (known as presence/absence rasters). In this case rasterize() requires only one argument in addition to x and y (the aforementioned vector and raster objects): a value to be transferred to all non-empty cells specified by field.

ch_raster1 = rasterize(cycle_hire_osm_projected, raster_template, field = 1)

The fun argument specifies summary statistics used to convert multiple observations in close proximity into associate cells in the raster object. By default fun = "last" is used but other options such as fun = "count" can be used, in this case to count the number of cycle hire points in each grid cell.

ch_raster2 = rasterize(cycle_hire_osm_projected, raster_template, 
                       field = 1, fun = "count")

The new output, ch_raster2, shows the number of cycle hire points in each grid cell. The cycle hire locations have different numbers of bicycles described by the capacity variable, raising the question, what’s the capacity in each grid cell? To calculate that we must sum the field ("capacity"), resulting in output illustrated in panel (D) below, calculated with the following command (other summary functions such as mean could be used):

ch_raster3 = rasterize(cycle_hire_osm_projected, raster_template, 
                       field = "capacity", fun = sum)
Examples of point rasterization.

Examples of point rasterization.

Another dataset based on California’s polygons and borders (created below) illustrates rasterization of lines. After casting the polygon objects into a multilinestring, a template raster is created with a resolution of a 0.5 degree:

california = dplyr::filter(us_states, NAME == "California")
california_borders = st_cast(california, "MULTILINESTRING")
raster_template2 = raster(extent(california), resolution = 0.5,
                         crs = st_crs(california)$proj4string)

Line rasterization is demonstrated in the code below. In the resulting raster, all cells that are touched by a line get a value.

california_raster1 = rasterize(california_borders, raster_template2)

Polygon rasterization, by contrast, selects only cells whose centroids are inside the selector polygon.

california_raster2 = rasterize(california, raster_template2)
Examples of line and polygon rasterizations.

Examples of line and polygon rasterizations.

As with raster::extract(), raster::rasterize() works well for most cases but is not performance optimized. Fortunately, there are several alternatives, including the fasterize::fasterize() and gdalUtils::gdal_rasterize(). The former is much (100 times+) faster than rasterize(), but is currently limited to polygon rasterization. The latter is part of GDAL and therefore requires a vector file (instead of an sf object) and rasterization parameters (instead of a Raster* template object) as inputs.1

Spatial vectorization

Spatial vectorization is the counterpart of rasterization, but in the opposite direction. It involves converting spatially continuous raster data into spatially discrete vector data such as points, lines or polygons.

Be careful with the wording! In R, vectorization refers to the possibility of replacing for-loops and alike by doing things like 1:10 / 2 (see also @wickham_advanced_2014).

The simplest form of vectorization is to convert the centroids of raster cells into points. rasterToPoints() does exactly this for all non-NA raster grid cells. Setting the spatial parameter to TRUE is needed to get a spatial object instead of a matrix.

elev_point = rasterToPoints(elev, spatial = TRUE) %>% 
  st_as_sf()
Raster and point representation of the elev object.

Raster and point representation of the elev object.

Another common type of vectorization involves conversion of rasters to polygons. This can be done with raster::rasterToPolygons(), which converts each raster cell into a polygon consisting of five coordinates, all of which are stored in memory (explaining why rasters are often fast compared with vectors!).

This is illustrated below by converting the grain object into polygons and subsequently dissolving borders between polygons with the same attribute values (also see the dissolve argument in rasterToPolygons()). Attributes in this case are stored in a column called layer. (Note: a convenient alternative for converting rasters into polygons is spex::polygonize() which by default returns an sf object.)

grain_poly = rasterToPolygons(grain) %>% 
  st_as_sf()
grain_poly2 = grain_poly %>% 
  group_by(layer) %>%
  summarize()
Illustration of vectorization of raster (left) into polygon (center) and polygon aggregation (right).

Illustration of vectorization of raster (left) into polygon (center) and polygon aggregation (right).

  1. See more at http://gdal.org/gdal_rasterize.html.↩